library(clusterProfiler)
simpl_enrich <- simplify
library(enrichplot)
library(tidyverse)
library(ggpubr)
library(patchwork)
library(ggrepel)
source('00_util_scripts/mod_bplot.R')

kegg_mva <-
  c('Acat1','Acat2','Hmgcs1','Hmgcs2','Hmgcr','Mvk','Pmvk','Mvd','Idi1','Idi2','Fdps')

mva_fct <- tibble(gene = kegg_mva, ordered = fct_inorder(kegg_mva))

kegg_d.mva <-
  c('FDFT1','LSS','SQLE','GGPS1','DOLK','PDSS1','PDSS2') |>
  str_to_title()

key_cytokine <- c('Il6','Ccl20','Tslp','Flt3l','Csf2','Tnf')

meta_kera <- read_csv('mission/FPP/zww_sa_mice/results/aureus_kera_meta.csv')

kera_infct_deg <- read_csv('mission/FPP/zww_sa_mice/results/kera.infect.deg.csv')

kera_v3_deg <-
  read_csv('mission/FPP/zww_sa_mice/results/mice_v3h_SA-PBS_deg.csv')

kera_v3l_deg <-
  read_csv('mission/FPP/zww_sa_mice/results/mice_v3l_SA-PBS_deg.csv')

sa_v3hvl_deg <-
  read_csv('mission/FPP/zww_sa_mice/results/mice_SA_v3h-v3l_deg.csv')

ct_v3hvl_deg <-
  read_csv('mission/FPP/zww_sa_mice/results/mice_ctrl_v3h-v3l_deg.csv')

# infect v3h vs ctrl v3h ------------
### FIG: SA v3h vs v3l kera volcano ===========
kera_v3_deg |>
  filter(!str_detect(gene, 'Rik$|^Gm|AC116')) |>
  plot_bill_volc('SA infection','PBS control',
                 highlights = c(kegg_mva, key_cytokine)) +
  ggtitle('Infected Trpv3-hi-KC vs control Trpv3-hi-KC')

kera_v3_deg |>
  plot_pub_volc('SA V3-hi KC','PBS V3-hi KC',
                 highlights = c(kegg_mva, key_cytokine))

publish_pdf('mission/FPP/micefig2/mice_v3h_SA-pbs_volc.pdf')

v3_inf_up <- kera_v3_deg |>
  filter(p_val_adj < .05 & avg_log2FC > 2)

v3_inf_down <- kera_v3_deg |>
  filter(p_val_adj < .05 & avg_log2FC < -2)

v3_inf_up_ora_gosim <- v3_inf_up$gene |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T) |>
  simpl_enrich() 

v3_inf_up_ora_gosim |>
  pluck('result') |>
  head(12) |>
  plot_enrichment() +
  ggtitle('Upregulated GO BP pathway in Trpv3+ keratinocytes after infection')

v3_inf_down_ora <- v3_inf_down$gene |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T) |>
  simpl_enrich() 

v3_inf_down_ora |>
  pluck('result') |>
  head(12) |>
  plot_enrichment('blue') +
  ggtitle('Downregulated GO BP pathway in Trpv3+ keratinocytes after infection')

v3_inf_up_ora_gosim@result |>
  as_tibble() |>
  filter(qvalue < .05) |>
  write_csv('mission/FPP/v3_inf-vs-pbs_up_go.csv')

v3_inf_down_ora@result |>
  as_tibble() |>
  filter(qvalue < .05) |>
  write_csv('mission/FPP/v3_inf-vs-pbs_down_go.csv')

### FIG: SA vs ctrl V3h-KC upgo ===========
v3_inf_up_ora_godf <-
  read_csv('mission/FPP/zww_sa_mice/results/v3_inf-vs-pbs_up_go.csv')

v3_inf_up_ora_godf |>
  publish_enrichment() +
  ggtitle('Upregulated GO BP pathway in V3-hi-KC: SA vs PBS')

publish_pdf('mission/FPP/micefig2/mice_v3h_SA-pbs_upgo.pdf', width = 60)

v3_inf_down_ora_godf <-
  read_csv('mission/FPP/v3_inf-vs-pbs_down_go.csv')

v3_inf_down_ora_godf |>
  head(12) |>
  plot_enrichment('blue') +
  ggtitle('Downregulated GO BP pathway in Trpv3+ keratinocytes after infection')

ggsave('mice_v3h_pbs-SA_downgo.pdf',
       width = 7, height = 4)

## interested genes in pathways? ---------
v3_inf_up_ora_godf |>
  separate_longer_delim(geneID, delim = '/') |>
  filter(geneID %in% chole_p_str) |>
  select(Description, qvalue)

v3_inf_up_ora_godf |>
  separate_longer_delim(geneID, delim = '/') |>
  filter(geneID %in% key_cytokine) |>
  select(Description, qvalue, geneID)

leuk_migr <- v3_inf_up_ora_godf |>
  separate_longer_delim(geneID, delim = '/') |>
  filter(Description == 'leukocyte migration') |>
  pull(geneID)

v3_inf_up |>
  filter(gene %in% leuk_migr)

sobj_kera |>
  filter(trpv3_status == 'Trpv3-high') |>
  ScaleData(features = leuk_migr) |>
  DoHeatmap(features = leuk_migr,
            group.by = 'orig.ident',
            angle = 0,
            vjust = 1,
            hjust = .5) +
  scale_fill_gradient2(high = 'red', low = 'blue')

g1 <- last_plot()

g1[[1]] + labs(title = 'GO BP pathway: Leukocyte migration')

ggsave('mice_v3h_SAup_GO1_heatmap.pdf',
       width = 6, height = 4)

small_catabo <- v3_inf_up_ora_godf |>
  separate_longer_delim(geneID, delim = '/') |>
  filter(Description == 'small molecule catabolic process') |>
  pull(geneID)

sobj_kera |>
  filter(trpv3_status == 'Trpv3-high') |>
  ScaleData(features = small_catabo) |>
  DoHeatmap(features = small_catabo,
            group.by = 'orig.ident',
            angle = 0,
            vjust = 1,
            hjust = .5) +
  scale_fill_gradient2(high = 'red', low = 'blue')

g2 <- last_plot()

g2[[1]] + labs(title = 'GO BP pathway: Small molecule catabolic process')

ggsave('mice_v3h_SAup_GO2_heatmap.pdf',
       width = 6, height = 4)

# infect v3l vs ctrl v3l ------------
### FIG: SA v3l vs PBS v3l kera volcano ===========
kera_v3l_deg |>
  filter(!str_detect(gene, 'Rik$|^Gm|AC116'), p_val_adj != 0) |>
  plot_bill_volc('SA V3-lo KC','PBS V3-lo KC',
                 highlights = c(kegg_mva, key_cytokine))

kera_v3l_deg |>
  filter(p_val_adj != 0) |>
  plot_pub_volc('SA V3-lo KC','PBS V3-lo KC',
                 highlights = c(kegg_mva, key_cytokine))

publish_pdf('mission/FPP/micefig2/mice_v3l_SA-pbs_volc.pdf')

v3l_inf_up <- kera_v3l_deg |>
  filter(p_val_adj < .05 & avg_log2FC > 1)

v3l_inf_down <- kera_v3l_deg |>
  filter(p_val_adj < .05 & avg_log2FC < -1)

v3l_inf_upgo <- v3l_inf_up$gene |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T) |>
  simpl_enrich() 

v3l_inf_upgo@result |>
  plot_enrichment() +
  ggtitle('Upregulated GO BP pathway in Trpv3-lo keratinocytes after infection')

v3l_inf_upgo@result |>
  publish_enrichment() +
  ggtitle('Upregulated GO BP pathway in V3-lo KC: SA vs PBS')

publish_pdf('mission/FPP/micefig2/v3l_sa_vs_pbs_upgo.pdf', width = 65)

v3_inf_down_ora <- v3_inf_down$gene |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T) |>
  simpl_enrich() 

v3_inf_down_ora |>
  pluck('result') |>
  head(12) |>
  plot_enrichment('blue') +
  ggtitle('Downregulated GO BP pathway in Trpv3+ keratinocytes after infection')

v3_inf_up_ora_gosim@result |>
  as_tibble() |>
  filter(qvalue < .05) |>
  write_csv('mission/FPP/v3_inf-vs-pbs_up_go.csv')

v3_inf_down_ora@result |>
  as_tibble() |>
  filter(qvalue < .05) |>
  write_csv('mission/FPP/v3_inf-vs-pbs_down_go.csv')

### FIG: SA vs ctrl V3h-KC upgo ===========
v3_inf_up_ora_godf <-
  read_csv('mission/FPP/zww_sa_mice/results/v3_inf-vs-pbs_up_go.csv')

v3_inf_up_ora_godf |>
  head(10) |>
  plot_enrichment() +
  ggtitle('Upregulated GO BP pathway in SA-infected vs control V3-hi-KC')

publish_pdf('micefig/mice_v3h_SA-pbs_upgo.pdf', width = 60)

v3_inf_down_ora_godf <-
  read_csv('mission/FPP/v3_inf-vs-pbs_down_go.csv')

v3_inf_down_ora_godf |>
  head(12) |>
  plot_enrichment('blue') +
  ggtitle('Downregulated GO BP pathway in Trpv3+ keratinocytes after infection')

# infect kera vs ctrl kera -------
### FIG: SA vs ctrl kera volcano ===========
kera_infct_deg |>
  filter(gene %in% c(key_cytokine, kegg_mva))

kera_infct_deg |>
  filter(!str_detect(gene, 'Rik$|^Gm')) |>
  plot_pub_volc('keratinocytes', highlights = c(key_cytokine, kegg_mva))

publish_pdf('micefig/mice_pbs-kera_SA-kera_volcano.pdf')

### FIG: SA vs ctrl kera upgo ===========
kr_inf_up <- kera_infct_deg |>
  filter(p_val_adj < .05 & avg_log2FC > 2)

kr_inf_down <- kera_infct_deg |>
  filter(p_val_adj < .05 & avg_log2FC < -2)

kr_inf_up_ora_go <- kr_inf_up$gene |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T) |>
  simpl_enrich() 

kr_inf_up_ora_go |>
  pluck('result') |>
  head(12) |>
  plot_enrichment() +
  ggtitle('Upregulated GO BP pathway in all keratinocytes after infection')

kr_inf_down_ora <- kr_inf_down$gene |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T) |>
  simpl_enrich() 

kr_inf_down_ora |>
  pluck('result') |>
  head(12) |>
  plot_enrichment('blue') +
  ggtitle('Downregulated GO BP pathway in all keratinocytes after infection')

kr_inf_up_ora_go@result |>
  as_tibble() |>
  filter(qvalue < .05) |>
  write_csv('mission/FPP/kr_inf-vs-pbs_up_go.csv')

kr_inf_down_ora@result |>
  as_tibble() |>
  filter(qvalue < .05) |>
  write_csv('mission/FPP/kr_inf-vs-pbs_down_go.csv')

read_csv('mission/FPP/kr_inf-vs-pbs_up_go.csv') |>
  head(10) |>
  publish_enrichment()

publish_pdf('micefig/mice_all_kera_pbs-SA_upgo.pdf', width = 65)

read_csv('mission/FPP/kr_inf-vs-pbs_down_go.csv') |>
  head(12) |>
  plot_enrichment('blue') +
  ggtitle('Downregulated GO BP pathway in all keratinocytes after infection')

ggsave('mice_all_kera_pbs-SA_downgo.pdf',
       width = 7, height = 4)

# compare v3h vs v3l ---------
## FIG: SA v3h vs v3l volcano =============
sa_v3hvl_deg |>
  filter(!str_detect(gene, 'AC116') & p_val_adj != 0) |>
  plot_bill_volc('V3-hi KC',group2 = 'V3-lo KC',
                 highlights = c(key_cytokine, kegg_mva))

sa_v3hvl_deg |>
  filter(!str_detect(gene, 'AC116') & p_val_adj != 0) |>
  plot_pub_volc('V3-hi KC',group2 = 'V3-lo KC',
                 highlights = c(key_cytokine, 'Mvk'))

publish_pdf('mission/FPP/micefig2/mice_SA_v3h-v3l_volcano.pdf')

inf_v3hvl_cyto <- trpv3_kera_marker |>
  filter(gene %in% key_cytokine)

inf_v3hvl_cyto |>
  ggplot(aes(x = 'Trpv3-hi vs Trpv3-lo', y = gene, color = avg_log2FC, size = -log10(p_val_adj))) +
  geom_point() +
  scale_color_gradient2(high = 'red',low = 'blue') +
  theme_pubr(legend = 'right') +
  labs_pubr() +
  labs(x = '')

ggsave('mice_SA_v3h-v3l-logfc_cytokine_bubbleplot.pdf',
       width = 4, height = 3.5)

## GO ORA ---------
v3_mrkr_up <- sa_v3hvl_deg |>
  filter(p_val_adj < .05 & avg_log2FC > 1)

v3_mrkr_down <- sa_v3hvl_deg |>
  filter(p_val_adj < .05 & avg_log2FC < -1)

v3_mrkr_up_go <- v3_mrkr_up$gene |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T) |>
  simpl_enrich() 

v3_mrkr_down_go <- v3_mrkr_down$gene |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T) |>
  simpl_enrich() 

v3_mrkr_up_go@result |>
  publish_enrichment() +
  ggtitle('Upregulated GO BP pathway in SA group: V3-hi KC vs V3-lo KC')

publish_pdf('mission/FPP/micefig2/mice_SA_v3hvl_upgo.pdf', width = 65)

v3_mrkr_down_go@result |>
  head(12) |>
  plot_enrichment('blue') +
  ggtitle('Downregulated GO BP pathway in Trpv3+ vs Trpv3- keratinocytes')

v3_mrkr_up_go@result |>
  as_tibble() |>
  filter(qvalue < .05) |>
  write_csv('mission/FPP/inf_v3h-vs-v3l_up_go.csv')

v3_mrkr_down_go@result |>
  as_tibble() |>
  filter(qvalue < .05) |>
  write_csv('mission/FPP/inf_v3h-vs-v3l_down_go.csv')

## FIG: pbs v3h vs v3l volcano =============
ct_v3hvl_deg |>
  filter(!str_detect(gene, 'AC116') & p_val_adj != 0) |>
  plot_bill_volc('PBS Trpv3-hi KC',group2 = 'PBS Trpv3-lo KC',
                 highlights = c(key_cytokine,kegg_mva)) +
  ggtitle('PBS Trpv3-hi KC vs PBS Trpv3-lo KC')

ct_v3hvl_deg |>
  filter(!str_detect(gene, 'AC116') & p_val_adj != 0) |>
  plot_pub_volc('PBS V3-hi KC',group2 = 'PBS V3-lo KC',
                 highlights = c(key_cytokine,kegg_mva))

publish_pdf('mission/FPP/micefig2/mice_pbs_v3h-v3l_volcano.pdf')

## GO ORA ---------
ct_v3hvl_up <- ct_v3hvl_deg |>
  filter(p_val_adj < .05 & avg_log2FC > 1)

ct_v3hvl_down <- ct_v3hvl_deg |>
  filter(p_val_adj < .05 & avg_log2FC < -1)

ct_v3hvl_upgo <- ct_v3hvl_up$gene |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T) |>
  simpl_enrich() 

v3_mrkr_down_go <- v3_mrkr_down$gene |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T) |>
  simpl_enrich() 

ct_v3hvl_upgo@result |>
  plot_enrichment() +
  ggtitle('Upregulated GO BP pathway in control group Trpv3-hi vs Trpv3-lo KC')

ct_v3hvl_upgo@result |>
  publish_enrichment() +
  ggtitle('Upregulated GO BP pathway in PBS group: V3-hi KC vs V3-lo KC')

publish_pdf('mission/FPP/micefig2/pbs_v3hvl_upgo.pdf', width = 65)

ct_v3hvl_upgo |>
  as_tibble() |>
  filter(qvalue < .05) |>
  write_csv('mission/FPP/zww_sa_mice/pbs_v3hvl_upgo.csv')

read_tsv('mission/FPP/zww_sa_mice/results/pbs.v3hvl.upgo.digest.tsv') |>
  plot_enrichment() +
  ggtitle('Upregulated GO pathway in PBS group v3h vs v3l KC')

v3_mrkr_down_go@result |>
  as_tibble() |>
  filter(qvalue < .05) |>
  write_csv('mission/FPP/inf_v3h-vs-v3l_down_go.csv')

### FIG: SA v3h vs v3l upgo ============
read_csv('mission/FPP/inf_v3h-vs-v3l_up_go.csv') |>
  head(10) |>
  publish_enrichment()

publish_pdf('micefig/mice_SA_kera_v3h-v3l_upgo.pdf', width = 60)

### FIG: v3l SA-PBS upgo ===============
v3l_sa_upgo <- read_csv('mission/FPP/v3l_kera_SA-PBS_deg.csv') |>
  filter(p_val_adj < .05 & avg_log2FC > 2) |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           readable = T) |>
  simpl_enrich() 

v3l_sa_upgo@result |>
  head(10) |>
  publish_enrichment()

publish_pdf('micefig/mice_SA_v3l_SA-PBS_upgo.pdf', width = 60)

read_csv('mission/FPP/inf_v3h-vs-v3l_down_go.csv') |>
  head(12) |>
  plot_enrichment('blue') +
  ggtitle('Downregulated GO BP pathway in Trpv3+ vs Trpv3- keratinocytes')

ggsave('mice_SA_kera_v3h-v3l_downgo.pdf',
       width = 7, height = 4)

v3_mrkr_up_go@result |>
  as_tibble() |>
  filter(str_detect(Description, 'retinoid meta')) |>
  pull(geneID)

# compare cell frac change upon infection --------
stat_kera <- meta_kera |>
  count(orig.ident, seurat_clusters) |>
  group_by(orig.ident) |>
  calc_frac_conf_on_grouped_count()

pval_kera <- meta_kera |>
  test_on_grouped_count(orig.ident, seurat_clusters) |>
  ungroup() |>
  mutate(p.adj = p.adjust(p.value, method = 'BH'),
         seurat_clusters = subtype)

stat_kera |>
  ungroup() |>
  select(orig.ident, seurat_clusters, fraction) |>
  pivot_wider(names_from = orig.ident, values_from = fraction) |>
  left_join(pval_kera) |>
  mutate(fraction_log2fc = log2(infected / PBS),
         type = case_when(p.adj < .05 & fraction_log2fc < 0 ~ 'Decreased after infection',
                          p.adj < .05 & fraction_log2fc > 0 ~ 'Increased after infection',
                          .default = 'NS'),
         clusters = as.character(seurat_clusters) |> fct_reorder(seurat_clusters)) |>
  ggplot(aes(clusters, fraction_log2fc, fill = type)) +
  geom_col() +
  coord_flip() +
  theme_pubr(legend = 'right') +
  scale_fill_manual(values = c('blue','red','grey')) +
  labs(y = 'Log2 fold changes after infection',
       title = 'Cell fraction changes of keratinocytes after infection')

## FIG: v3h frac in PBS & SA ===========
v3_frac_conf <- meta_kera |>
  calc_frac_conf_on_grouped_count(orig.ident, trpv3_status)

v3_frac_pval <- meta_kera |>
  test_on_grouped_count(orig.ident, trpv3_status)

v3_frac_conf |>
  filter(trpv3_status == 'Trpv3-high') |>
  mutate(group = fct_relevel(orig.ident, 'PBS')) |>
  ggplot(aes(group, fraction, ymin = conf.low, ymax = conf.high, fill = group)) +
  geom_col() +
  geom_errorbar(width = .5) +
  scale_fill_manual(values = c('blue','red'), labels = c('PBS','SA')) +
  labs(y = 'Fraction in all keratinocytes', x = 'Group',
       title = 'Fraction of Trpv3-hi KC in mice skin') +
  theme_pubr(legend = 'right') +
  theme_jpub

publish_pdf('mission/FPP/micefig2/mice_v3h_krea_SA-PBS_fraction_barplot.pdf',
            width = 40)

# compare 8 vs 10 kera ----------
kera8v10 <- read_csv('mission/FPP/zww_sa_mice/kera8v10.deg.csv')

key_signaling <- c('Ptk2b','Kras','Nfkb1','Nfkb2','Calm1','Ppp3cc','Nfam1')

kera8v10 |>
  filter(gene %in% key_signaling)

## go ora ------------
kera8_upgo <- kera8v10 |>
  filter(avg_log2FC > 2 & p_val_adj < .05) |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           qvalueCutoff = .05,
           readable = T)
  
kera8_upgo@result |>
  head(12) |>
  plot_enrichment() +
  ggtitle('Upregulated GO pathways in cluster 8 keratinocytes')

kera8_upgo@result |>
  as_tibble() |>
  filter(p.adjust < .05) |>
  DT::datatable()

kera10_upgo <- kera8v10 |>
  filter(avg_log2FC < 2 & p_val_adj < .05) |>
  pull(gene) |>
  enrichGO(OrgDb = 'org.Mm.eg.db',
           keyType = 'SYMBOL',
           ont = 'BP',
           qvalueCutoff = .05,
           readable = T)

kera10_upgo@result |>
  head(12) |>
  plot_enrichment() +
  ggtitle('Upregulated GO pathways in cluster 10 keratinocytes')

## go gsea ---------------
sig8v10 <- kera8v10 |>
  filter(p_val_adj < .05)

gsgo8v10 <- sig8v10$avg_log2FC |>
  set_names(sig8v10$gene) |>
  sort(decreasing = T) |>
  gseGO(OrgDb = 'org.Mm.eg.db',
        keyType = 'SYMBOL',
        ont = 'BP',
        pvalueCutoff = 1)

gsgo8v10@result |>
  as_tibble() |>
  filter(str_detect(Description, 'alcium')) |>
  DT::datatable()

gsgo8v10 |>
  gseaplot2(geneSetID = c('GO:0071277','GO:0019722'), pvalue_table = T)

gsgo10v8 <- -sig8v10$avg_log2FC |>
  set_names(sig8v10$gene) |>
  sort(decreasing = T) |>
  gseGO(OrgDb = 'org.Mm.eg.db',
        keyType = 'SYMBOL',
        ont = 'BP')

gsgo10v8@result |>
  as_tibble() |>
  slice_max(NES, n =3)

# pathview -----
v3h.sa.vs.pbs <- kera_v3_deg |>
  pull(avg_log2FC) |>
  set_names(kera_v3_deg$gene)

kera_v3_deg |>
  filter(gene %in% kegg_d.mva) |>
  mutate(gene = fct_relevel(gene, 'Fdft1', 'Sqle','Lss') |> fct_rev()) |>
  ggplot(aes('', gene, fill = avg_log2FC)) +
  geom_tile() +
  scale_fill_gradient2(low = 'blue', high = 'red', limits = c(-1,2),
                       na.value = 'grey80') +
  theme_pubr(legend = 'right') +
  labs(x = '', y = 'Gene')

kera_v3_deg |>
  right_join(mva_fct) |>
  ggplot(aes('', fct_rev(ordered), fill = avg_log2FC)) +
  geom_tile() +
  scale_fill_gradient2(low = 'blue', high = 'red', limits = c(-1,NA),
                       na.value = 'grey80') +
  theme_pubr(legend = 'right') +
  labs(x = '', y = 'Gene')

library(pathview)
png.out <- pathview(
  gene.data = v3h.sa.vs.pbs,
  gene.idtype = 'SYMBOL',
  species = 'mouse',
  pathway.id = '00900',
  out.suffix = "v3h.sa.vs.pbs",
  kegg.native = T
)

pdf.out <- pathview(
  gene.data = v3h.sa.vs.pbs,
  gene.idtype = 'SYMBOL',
  species = 'mouse',
  pathway.id = '00900',
  out.suffix = "v3h.sa.vs.pbs",
  kegg.native = F,
  rankdir = 'TB',
  same.layer = F,
  pdf.size = c(14,7),
  expand.node = T
)

pv.out2 <- keggview.graph(
  gene.data = v3h.sa.vs.pbs,
  gene.idtype = 'SYMBOL',
  species = 'mouse',
  pathway.id = '00900',
  out.suffix = "v3h.sa.vs.pbs",
  kegg.native = F,
  plot.data.gene = kegg_mva,
  path.graph = pv.out$plot.data.gene
)
